home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / FLTRBTTR.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  4KB  |  116 lines
  1. /* 
  2.  * fltrbttr.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. #include <stdlib.h>
  10. #include <malloc.h>
  11. #include <math.h>
  12.  
  13. /*
  14.  * fltrbttr()
  15.  *   DESCRIPTION:
  16.  *     fltrbttr performs Butterworth filtering in frequency domain
  17.  *     with specified low & high pass bounds and filter order
  18.  *   ARGUMENTS:
  19.  *     imgReal(float **) pointer to real image array
  20.  *     imgImag(float **) pointer to imaginary image array
  21.  *     nRow, nCol(long) number of rows and columns for real and img arrays
  22.  *     passType(short)  passType=0 low pass
  23.  *                              passType=1 high pass
  24.  *                              passType=2 band pass
  25.  *                              passType=3 stop pass
  26.  *     f1, f2(double) low, high pass frequency bounds
  27.  *     order(double) order of filter
  28.  *   RETURN VALUE:
  29.  *     none
  30.  */
  31.  
  32. void
  33. fltrbttr (imgReal, imgImag, nRow, nCol, passType, f1, f2, order)
  34.      float **imgReal, **imgImag;  /* real,complex image arrays */
  35.      long nRow, nCol;           /* size of image */
  36.      short passType;            /* lowpass=0, high=1, band=2, stop=3 */
  37.      double f1, f2;             /* low,high pass frequency bounds */
  38.      double order;              /* order of filter */
  39. {
  40.   register long x, y, nRowD2, nColD2,  /* half of row/column */
  41.     nRowM1, nColM1;             /* row/col minus 1 */
  42.   double bigger;                /* passband fraction fct of bigger axis */
  43.   double *fltr;                 /* filter vector */
  44.   long maxRad;                  /* max. sampling radius */
  45.   long dist;                    /* euclidean distance to point */
  46.   double tempL, tempH;
  47.   long i;
  48.  
  49. /* allocate filter vector -- radius from (u,v)=(0,0) to corners of freq plot */
  50.   nRowD2 = nRow / 2;
  51.   nColD2 = nCol / 2;
  52.   bigger = (double) (nRowD2 > nColD2) ? nRowD2 : nColD2;
  53.   maxRad = (long) sqrt ((double) (nRowD2 * nRowD2 + nColD2 * nColD2)) + 1;
  54.   fltr = (double *) calloc (maxRad, sizeof (double));
  55.  
  56. /* determine filter vector coefficients for chosen passband type of filter */
  57.   order *= 2.0;
  58.   switch (passType) {
  59.   case 0:                      /* low-pass */
  60.     f1 *= bigger;
  61.     for (i = 0; i < maxRad; i++) {
  62.       fltr[i] = 1.0 / (1.0 + 0.414 * pow (((double) i / f1), order));
  63.  
  64.     }
  65.     break;
  66.   case 1:                      /* high-pass */
  67.     f1 *= bigger;
  68.     for (i = 0; i < maxRad; i++)
  69.       fltr[i] = 1.0 - 1.0 / (1.0 + 0.414 * pow (((double) i / f1), order));
  70.     break;
  71.   case 2:                      /* band-pass */
  72.     f1 *= bigger;
  73.     f2 *= bigger;
  74.     for (i = 0; i < maxRad; i++) {
  75.       tempL = 1.0 / (1.0 + 0.414 * pow (((double) i / f1), order));
  76.       tempH = 1.0 - 1.0 / (1.0 + 0.414 * pow (((double) i / f2), order));
  77.       fltr[i] = tempL * tempH;
  78.     }
  79.     break;
  80.   case 3:                      /* and band-stop */
  81.     f1 *= bigger;
  82.     f2 *= bigger;
  83.     for (i = 0; i < maxRad; i++) {
  84.       tempL = 1.0 / (1.0 + 0.414 * pow (((double) i / f1), order));
  85.       tempH = 1.0 - 1.0 / (1.0 + 0.414 * pow (((double) i / f2), order));
  86.       fltr[i] = tempL + tempH;
  87.     }
  88.     break;
  89.   default:
  90.     break;
  91.   }
  92.  
  93. /* perform filtering */
  94.   nRowM1 = nRow - 1;
  95.   nColM1 = nCol - 1;
  96.   for (y = 0; y < nRowD2; y++) {
  97.     for (x = 0; x < nColD2; x++) {
  98.       dist = (long) (sqrt ((double) (y * y + x * x)) + 0.5);
  99.  
  100.  
  101.       imgReal[y][x] *= (float) fltr[dist];
  102.       imgReal[y][nColM1 - x] *= (float) fltr[dist];
  103.       imgReal[nRowM1 - y][x] *= (float) fltr[dist];
  104.       imgReal[nRowM1 - y][nColM1 - x] *= (float) fltr[dist];
  105.  
  106.       imgImag[y][x] *= (float) fltr[dist];
  107.       imgImag[y][nColM1 - x] *= (float) fltr[dist];
  108.       imgImag[nRowM1 - y][x] *= (float) fltr[dist];
  109.       imgImag[nRowM1 - y][nColM1 - x] *= (float) fltr[dist];
  110.     }
  111.   }
  112.  
  113.  
  114.   return;
  115. }
  116.